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RESUMEN 

El comportamiento de tanques para almacenamiento de Ifquidos peligrosos ha 
demostrado su vulnerabilidad frente a terremotos destructives. Su colapso conduce 
al derrame de las sustancias contenidas (inflamables, explosivas, toxicas) 
impactando negativamente en el medio ambiente circundante y provocando grandes 
perdidas economicas para la region donde se encuentran emplazados. El objetivo 
de este trabajo es determinar la efectividad de los sistemas de aislacion sfsmica en 
tanques cilfnd ricos apoyados que almacenan Ifquidos peligrosos, y medir la 
incidencia de determinados parametros sobre la respuesta. Para lograr esto, es 
necesario realizar un analisis parametrico de la historia de la respuesta, con 
multiples corridas y un costo computacional elevado. Por ello, en este artfculo se 
propone un metodo simple y rapido, basado en la resolution del problema a traves 
de la Ecuacion de Euler-Lagrange. Como dispositivos de aislamiento se utilizaron 
aisladores friccionales de simple pendulo de friccion. Como input se utilizaron 
terremotos con caracterfsticas de movimientos sfsmicos impulsivos. Se modelaron 
un tanque con base fija y otro con base aislada, y se usaron como parametros de 
estudio la relation de esbeltez del tanque y, el perfodo y el coeficiente de friccion 
del sistema de aislacion sfsmica. Se determinaron las respuestas sfsmicas en 
terminos de desplazamientos de la superficie libre del Ifquido contenido, 
desplazamientos del sistema de aislacion y corte basal normalizado. Los resultados 
indican que el sistema de aislamiento reduce la respuesta maxima en todos los 
casos, pero la misma esta fuertemente condicionada por las caracterfsticas del 
registro impulsivo. 
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ABSTRACT 

The behavior of tanks for hazardous liquid storage has demonstrated their 
vulnerability to destructive earthquakes. Their collapse leads to the spillage of the 
contained substances (flammable, explosive, toxic) negatively impacting in the 
surrounding environment and causing big economic losses for the region where they 
are located. This work's goal is to determine the effectiveness of the systems of 
seismic isolation in supported cylindrical tanks, which store hazardous liquids, and 
to measure the incidence of certain parameters over the response. In order to 
achieve this, it is necessary to perform a parametric analysis of the response history, 
with multiple runs and with a high computational cost. Due to that, a simple and fast 
method, based on the resolution of the problem through an Euler-Lagrange Equation 
is proposed in this article. As isolation devices, friction isolators of single friction 
pendulum were used. As an input, earthquakes with impulsive seismic movement 
characteristics were used. Both a base fixed and a base isolated tank were modeled, 
and as study parameters, the slenderness ratio of the tank and the period and the 
friction coefficient of the seismic isolation system were modeled. The seismic 
responses in terms of displacements of the free surface of the contained liquid, 
displacements of the isolation system and the normalized base shear were 
determined. The results state that the isolation system reduces the maximum 
response in all cases, but it is strongly conditioned by the characteristics of the 
impulsive record. 

Keywords: seismic isolation, seismic response, cylindrical tank 
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1 INTRODUCCION 

El comportamiento de tanques para almacenamiento de Ifquidos peligrosos 
han demostrado su vulnerabilidad frente a terremotos destructives, tales como los 
ocurridos en Valdivia, Chile en 1960; Alaska, Estados Unidos en 1964; Northdrige, 
California, Estados Unidos, en 1994; Kobe, Japon en 1995; y Kocaeli, Turqufa en 
1999, entre otros. Su colapso conduce al derrame de las sustancias contenidas 
impactando negativamente en el medio ambiente circundante y provocando grandes 
perdidas economicas para la region donde se encuentran emplazados. Esto motiva 
el estudio del comportamiento sismico de los tanques, y el desarrollo de nuevas 
tecnicas para su protection, tales como la aislacion sfsmica de base. 

El comportamiento sismico de tanques ha sido estudiado desde hace varios 
anos por muchos autores. En 1933, Westergaard [1] determino la distribution de 
presiones de un fluido incompresible sobre una presa rfgida de pared vertical 
durante un terremoto. Posteriormente, en 1949 Jacobsen [2] determino la velocidad 
potencial de un tanque con fluido dentro y alrededor del mismo cuando experimenta 
en su base un desplazamiento traslacional impulsivo, despreciando la parte 
convectiva de la aceleracion. Anos mas tarde, en 1957 y 1963, Housner [3, 4] 
propuso ecuaciones simplificadas para contenedores de distintas formas y para ello 
realizo un analisis de la presion hidrodinamica desarrollada cuando un fluido 
contenido en un recipiente de paredes rfgidas es sujeto a aceleraciones 
horizontales, incluyendo presiones impulsivas y convectivas. En 1973, Veletsos [5] 
desarrollo un procedimiento simplificado para evaluar las fuerzas dinamicas 
inducidas por la componente lateral de un terremoto, teniendo en cuenta los efectos 
de flexibilidad del tanque, pero considerando solo las fuerzas impulsivas. Luego, en 
1981 y 1983, Haroun [6, 7] presento un metodo basado en la superposicion de los 
modos de vibracion lateral libre para evaluar la influencia de la flexibilidad de la 
pared en la respuesta sfsmica de un tanque. 

En cuanto a la tecnica del aislamiento sismico aplicada a tanques se puede 
citar los trabajos de Malhotra [8] los cuales demostraron que la respuesta sfsmica 
de tanques aislados se reduce considerablemente sin un incremento significativo 
del oleaje en la superficie libre. En 2001 , Wang [9] realizo simulaciones numericas 
de tanques con aisladores de pendulo de friccion (FPS), comprobando que su 
efectividad se incrementa conforme aumenta la relacion de esbeltez del tanque. En 
2002, Shrimali y Jangid [10] realizaron un analisis parametrico sobre un tanque 
ancho y otro esbelto, para medir la incidencia del perfodo de aislacion, el 
amortiguamiento y la fuerza de fluencia del sistema de aislacion sobre la respuesta 
pico del tanque, bajo una excitacion bidireccional considerando los efectos de 
interaccion, los cuales resultaron despreciables. En 2008 y 2012, Panchal y Jangid 
[11, 12] estudiaron la respuesta frente a sismos de falla cercana, de tanques 
aislados con dos dispositivos FPS novedosos, uno con coeficiente de friccion 
variable con el desplazamiento (VFPS) y otro con radio de curvatura variable con el 
desplazamiento (VCFPS), demostrando mejor comportamiento que los FPS con 
coeficiente de friccion y radio de curvatura constantes con el desplazamiento. En 
2010, Abali y Ugkan [13] realizaron un analisis parametrico de tanques con 
aisladores FPS, para medir la incidencia del perfodo de aislacion, la relacion de 
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aspecto del tanque y el coeficiente de friccion, en la respuesta del tanque, 
considerando la variacion de la carga axial sobre el sistema de aislacion debido al 
momento de vuelco y la aceleracion vertical, y encontrando que no debe ignorarse 
dicha variacion fundamentalmente en tanques esbeltos sometidos a sismos de falla 
cercana. En 2014, Saha [14] realizo un analisis parametrico para evaluar la 
respuesta de un tanque ancho y otro esbelto, utilizando aisladores modelados en 
forma lineal equivalente y en forma bilineal, mostrando que el modelo lineal 
equivalente sobreestima la respuesta de tanques esbeltos. 

El objetivo de este trabajo es evaluar la efectividad de los sistemas de 
aislacion sfsmica en tanques, y medir la incidencia de determinados parametros 
(periodo de aislacion, coeficiente de friccion y relacion de esbeltez) sobre la 
respuesta. Para lograr esto, es necesario realizar un analisis parametrico de la 
historia de la respuesta, con multiples corridas y un costo computacional elevado. 
Por ello, en este artfculo se propone un metodo simple y rapido, basado en la 
resolucion del problema a traves de la Ecuacion de Euler-Lagrange. 


2 METODOLOGIA 

Para lograr el objetivo planteado se siguen los siguientes pasos: i) elegir un 
modelo que represente el comportamiento ffsico del sistema: se usa el Modelo de 
Housner; ii) encontrar las ecuaciones diferenciales que gobiernan el sistema: se 
plantea la Ecuacion de Euler-Lagrange, que conduce a dos ecuaciones diferenciales 
ordinarias de segundo orden; iii) reducir el orden de las ecuaciones diferenciales: se 
usa la Ecuacion de Estado, que permite pasar de dos ecuaciones diferenciales 
ordinarias de segundo orden a cuatro ecuaciones diferenciales ordinarias de primer 
orden; iv) integrar las ecuaciones diferenciales: se emplea el metodo de integration 
numerica Runge Kutta de orden tres a paso constante. 

2.1 Modelo de Housner 

El modelo de Housner propone que la respuesta hidrodinamica de un 
sistema tanque-liquido se caracteriza por la participacion de dos contribuciones 
diferentes, llamadas componente impulsiva y componente convectiva. La 
componente impulsiva representa la porcion de Ifquido que se mueve al unfsono 
con las paredes del tanque. La componente convectiva representa el Ifquido que se 
mueve con un movimiento de oleaje de periodo grande en la porcion superior del 
tanque. Estas dos componentes pueden considerarse desacopladas, porque hay 
diferencias significativas en sus perfodos naturales. 

En un tanque apoyado sobre el suelo la amplitud del oleaje es un indicador 
de la intensidad del movimiento del suelo. Si un tanque con una superficie libre de 
agua se somete a una aceleracion horizontal del suelo X(t), las fuerzas ejercidas 
sobre el tanque por el agua son de dos tipos. Primero, cuando las paredes del 
tanque se aceleran, van y vuelven, y una parte del agua es forzada a participar en 
este movimiento, el cual ejerce una fuerza reactiva sobre el tanque, la misma que 
ejercerfa una masa mi conectada rfgidamente al tanque a una altura apropiada hi. 
La masa mi es fijada a la altura hi de modo que la fuerza horizontal ejercida por esta 
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sea colineal con la fuerza resultante ejercida por el agua equivalente. Segundo, el 
movimiento de la pared del tanque excita el agua en oscilaciones las cuales 
sucesivamente ejercen una fuerza oscilatoria sobre el tanque. Esta fuerza 
oscilatoria es la misma que ejercerfa una masa m c que puede oscilar 
horizontalmente contra un resorte de restriction. La masa m c corresponde al modo 
de oscilacion fundamental del agua, que es el modo mas importante en la mayorfa 
de los problemas sfsmicos. Si el sistema equivalente es sujeto a aceleraciones 
sfsmicas del suelo X(t), las fuerzas ejercidas sobre el tanque por las masas mi y m c 
seran las mismas que ejercerfa el agua sobre el tanque. En la figura 1 se 
esquematiza el modelo descripto. 



X(t) 

Figura 1 Modelo de Housner 

Para un tanque cilfndrico de radio R, y altura de Ifquido H, las ecuaciones 
para obtener los parametros necesarios del modelo son: 



tanhf 1 ^ 2 *) 

m = M v H J 

m 1 (1.732 R\ 

\ H ) 

( 1 ) 

m c 

*M 1S3 r h ) 

■ = 0 ' 835M (1.835 H) 

( 2 ) 
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( 3 ) 
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_ m c 2 g H 

oo c 2 m c = 4.032605 -f- 

c c MR 2 
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3 

h ^s H 

( 5 ) 

h c = H 

( R \ /0.9175 H\ I 

(l.835«) '“"'•I R )l 

( 6 ) 


Donde: R = radio del cilindro del tanque; H = altura del Ifquido; p = densidad 
del Ifquido; m = p n R 2 H masa del Ifquido contenido; mi = masa impulsiva; m c = 
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masa convectiva; g = 9.81 m/s 2 aceleracion de la gravedad; k c = rigidez convectiva; 
hi = altura impulsiva; h c = altura convectiva 

2.2 Ecuacion de Euler - Lagrange 

La Ecuacion de Euler - Lagrange se obtiene aplicando el Principio de los 
Trabajos Virtuales a la Segunda Ley de Newton. El Principio de Trabajos Virtuales 
(PTV) establece que “En un sistema mecanico, es condition necesaria y suficiente 
de equilibrio que el trabajo del conjunto de fuerzas aplicadas sobre desplazamientos 
virtuales compatibles con los vinculos sea nulo”. Para un sistema de M puntos 
materiales, el trabajo virtual total es: 

M 

8W = Y j Fi-$ r i = 0 (7) 

i-1 

Por otro lado, la segunda ley de Newton establece que: 

Ft = Pi (8) 

donde p t es el momentum lineal de la particula i-esima. De la expresion (8) 
se puede escribir: 


<Pi=Fi-Pi = 0 


(9) 


donde (p t se puede interpretar como la fuerza neta que habria que aplicar 
sobre la particula para mantenerla en “equilibrio”. Entonces, como en todo instante 
el sistema se mantiene en equilibrio bajo la accion de se cumple el principio de 
los trabajos virtuales: 


MM MM 

Fr Sri-^pr Sr t = 0 

i-1 i-1 i= 1 i=l 

Desarrollando el primer termino de (10) obtenemos: 



( 10 ) 


Z 

i = i 


Pi ■ = 




( 11 ) 


= ^ Qj dqj 

]= i 

donde q [Nxl] es el vector de grados de libertad del sistema (GDL), siendo 
N el numero de grados de libertad; y Qj = Z"i ( Pi-j , son las componentes de las 
fuerzas generalizadas. 


A continuacion, se desarrolla el segundo termino de (10): 
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Z pi ■ Sr> = Z mifi Z w, dq > = Z Z { mif ‘ ■ w) e< “ " 2) 

i= 1 i= 1 7=1 7 ;=1 i=l X J/ 


De la regia de la derivacion del producto interior, tenemos: 

da d db 

— ' b = —(a- b) — a- — 
dt dt dt 

Usando a = 7 x 1 , 1 -: y b = se obtiene: 

dqj 


NM 

ZZ{^ 


m,r 


i 1 1 


dr _ i 

dqj 


cL (dr: 


- m,r, — 


dt \dqj 


7=1 i = 1 

Si consideramos las siguientes igualdades: 


dqj 


d / dr i 


df t 

dqj 


dt \dqj } 

dr t dr t 
dqj dqj 

Entonces, reemplazando (15) y (16) en (14): 


N M 

w d 


dirt 
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— miTi 


dVi 
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(13) 


(14) 

(15) 

(16) 

(17) 


Ademas: 
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II 
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Donde T t = -mivf es la energia cinetica de la particula i-esima. Luego, 


reemplazando (18) y (19) en (17): 


V V(— (— \ 1 

n&Mydqj) dqj] 


dqj 


dT 


N , . 

:Y{±(— V 

4j [dt \dqjj dqj 


dqj 


(20) 


Reemplazando (11) en el primer termino de (10), y (20) en el segundo 
termino de (10), se obtiene: 
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( 22 ) 


En la ecuacion (22) Qj puede descomponerse como la suma de las fuerzas 
conservativas (derivadas de una funcion potencial) y las no conservativas: 


M 

dV dv; 

Reemplazando (23) en (22): 


(23) 


Zj| dt\dqj) dqj dqj 1 


Sqj = 0 


(24) 


Como los desplazamientos virtuales son arbitrarios, para que se cumpla la 
ecuacion (24) cada termino debe ser nulo, es decir: 


d ( dT\ dT dV 
dt \dqj ) dqj dqj 


+ - Qj c = 0 


(25) 


La ecuacion anterior representa una suma de fuerzas en equilibrio. Para 
aplicarla, se calcula la energia total (cinetica T y potencial V) del sistema (escalar) 
y se deriva respecto a las coordenadas generalizadas q (numero minimo de 
coordenadas independientes necesarias para describir el sistema = numero de 
grados de libertad), y se suman tambien las fuerzas no conservativas Q NC , para 
obtener de este modo las ecuaciones diferenciales que gobiernan el sistema. 

Los pasos a seguir para formular el problema en terminos de la Ecuacion 
de Euler - Lagrange son los siguientes: i) elegir un origen de coordenadas; ii) definir 
los grados de libertad; iii) definir vectores de position y velocidad; iv) calcular la 
energia cinetica; v) calcular la energia potencial (elastica mas gravitatoria); vi) 
calcular las fuerzas no conservativas; vii) armar la ecuacion. 

2.2.1 Origen de coordenadas 


El problema de un tanque con base fija y otro con base aislada, puede ser 
idealizado mediante los siguientes esquemas, a traves de la colocation de dos 
masas puntuales (impulsiva y convectiva) de acuerdo al Modelo de Housner. Los 
aisladores FPS pueden ser caracterizados por el radio de curvatura de una 
superficie concava, que define el periodo de aislacion, y por el coeficiente de friction 
en la superficie de deslizamiento, que define la disipacion de energia y 
amortiguamiento del sistema, y que si bien es variable en el tiempo con la velocidad 
de deslizamiento y la presion de contacto, para este caso se supone constante. El 
aislador se modela mediante una barra rigida con el mismo radio que la superficie 
concava de deslizamiento, dicha barra esta articulada en ambos extremos y 
conectada a otra barra horizontal rigida que representa la base del tanque, donde 
luego se empotra una barra vertical rigida que representa las paredes del tanque 
anclado, y que se conecta a las dos masas de Housner que representan el liquido 
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contenido. El origen de coordenadas se coloca en la union entre la barra vertical 
rigida (pared) y la barra horizontal rigida (fondo), tal como se muestra en la figura 2. 




Figura 2 Esquema simplificado del tanque: base fija (izquierda) y base aislada 

(derecha) 


2.2.2 Grados de libertad 

El problema del tanque aislado presenta dos grados de libertad, que se 
corresponden con los desplazamientos horizontales de las dos masas, la impulsiva 
(porcion de agua inferior y media que se mueve junto con las paredes rfgidas del 
tanque) con desplazamiento Xi generado por el movimiento pendular del sistema de 
aislacion al ser excitado por un sismo, y la convectiva (porcion de agua superior que 
oscila y genera oleaje en la superficie libre) con desplazamiento x c . Para el caso del 
tanque con base fija solo existe un grado de libertad o desplazamiento 
independiente, que corresponde a la masa convectiva, ya que la masa impulsiva al 
estar rfgidamente fija a las paredes del tanque y estas a su vez ancladas a la base, 
tiene un desplazamiento coincidente con el del suelo durante un terremoto. Esto se 
observa con claridad al dibujar los esquemas anteriores en posicion deformada, en 
la figura 3. 



Figura 3 Grados de libertad: base fija (izquierda) y base aislada (derecha) 

A partir de aquf se seguira el desarrollo solo para el caso del tanque aislado 
con dos grados de libertad. El caso del tanque con base fija y un solo grado de 
libertad puede obtenerse facilmente a partir del primero o incluso puede forzarse el 
modelo de base aislada para que funcione como base fija, colocando un radio de 
curvatura R c =0 y un coeficiente de friccion p=PGA. 
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Los grados de libertad, en terminos de desplazamiento, velocidad y 
aceleracion, son: 



(26) 

(27) 

(28) 


La rotacion 0 y el desplazamiento vertical Ay que generan los aisladores 
FPS debido a su movimiento pendular, pueden ser definidos en funcion del 
desplazamiento xi de la masa impulsiva. 



Figura 4 Desplazamientos y rotaciones en el aislador FPS 
A partir de la ecuacion de una circunferencia se despeja “y” : 

x? + y 2 = R 2 c (29) 


y = 



A partir de la figura 4 se obtienen Ay y 0: 


(30) 


Ay = R c -y = R c - R* - x\ 


send = — 
R r 


6 = arcsen 


© 


(31) 

(32) 

(33) 


2.2.3 Vectores de posicion y velocidad 


Se define la posicion de las masas impulsiva y convectiva, respecto al origen 
de coordenadas elegido, y en funcion de los grados de libertad asignados, llamando 
n al vector de posicion impulsivo y r c al convectivo. La coordenada “y” se obtiene 
como la altura h respectiva mas el desplazamiento vertical Ay que provocan los 
aisladores. 
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hi +R C - JR? -x? 


+ x c 


h c + R c - JR? - x? 


Los vectores de velocidad se obtienen a partir de los anteriores mediante el 
uso del Jacobiano J para obtener las derivadas de los vectores de posicion respecto 
a los grados de libertad. 

dr T 

f = —q=Jq ( 36 ) 

2.2.4 Energi'a cinetica 

Se calcula la energfa cinetica para cada masa del sistema con su respectivo 
vector de velocidad. 


T(q,q) =^^r£m fc r k =^^(J k q) T m k (J k q) 


2.2.5 Energi'a potencial 

Se calcula la energfa potencial gravitatoria V g (presente en las dos masas 
al levantarse por efecto de los aisladores) y la energfa potencial elastica V s (presente 
en el resorte de la masa convectiva). Luego se suman las dos contribuciones. 


Vg(q) = 'Y j rn k gr 2ik 


L, 


■u^m k u k 


2.2.6 Fuerzas no conservativas 

En este problema aparecen dos fuerzas no conservativas que disipan 
energfa, la fuerza F a de amortiguamiento en el movimiento convectivo del agua y la 
fuerza F r correspondiente al rozamiento en los aisladores de pendulo de friccion: 

P a = ~ c c v c sign(y c ) ( 40 ) 

c c = 2m c co c ^ c ( 41 ) 

F r = ~nNsign(Vi) ( 42 ) 

„ . mv t 2 \ „ . (vicose ) 2 ] 


N = mgcosQ H — — — = m + 
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Estas fuerzas se deben proyectar respecto a los grados de libertad, para 
ello se calcula el trabajo que realizan (vector position del punto de aplicacion por 
vector fuerza) y luego se deriva respecto a los grados de libertad. 

La fuerza de rozamiento del pendulo merece un analisis mas detallado. La 
fuerza normal interviniente varfa inversamente con el angulo de rotation 0 y es la 
suma de la componente radial del peso y la fuerza centrfpeta. El sentido de la fuerza 
de rozamiento es contrario al de la velocidad. 




l"m 


mg 


v t =0 




Figura 5 Magnitud de la fuerza de rozamiento y descomposicion de la velocidad 

en el piano (izquierda) 




Figura 6 Sentido de la fuerza de rozamiento y descomposicion de la fuerza en el 

piano (derecha) 
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2.2.7 Ecuacion de Euler Lagrange 

La Ecuacion de Euler-Lagrange se puede interpretar como la suma de 
distintos terminos de fuerzas en equilibrio, ya que las energfas del sistema se 
derivan respecto a los grados de libertad: 


d ( dT\ dT dV _ 
dt \dcjj) dqj + dqj 


A-B + C = D (44) 


d / dT\ _ d /dT\ .. d fdT\ . 
dt \dcjj) ~ dqj \dcjj) ^ + dqj \dcjj 1 

Agregando la fuerza sfsmica externa 


-> 4=41+42 (45) 

Fs a la suma de la fuerza total EL 


queda: 


EL = 41 +42 -B + C-D + Fs (46) 

En virtud que el sistema debe estar en equilibrio dinamico, la suma de 
todas las fuerzas internas y externas deber ser nula, por lo que: 

EL=A1+A2 — B + C — D + Fs = 0 (47) 

Asf se obtienen dos ecuaciones diferenciales de segundo orden, dado que 
el problema presenta dos grados de libertad. 

2.3 Ecuacion de Estado 


A partir de las dos ecuaciones de Euler-Lagrange obtenidas, se busca 
despejar la aceleracion de las masas impulsiva y convectiva. Para ello se arma una 
matriz lineal L colocando los terminos que acompanan a las aceleraciones 
buscadas. Esta matriz es la matriz de masa evaluada en q=0. Luego la parte no 
lineal de las ecuaciones se puede expresar como: 

NL = EL — a.L (48) 


Recordando que EL=0 se puede despejar la aceleracion como: 

a = — ZT 1 . NL (49) 



x c , %i> 

^2 (%ii X c , X„ X*c) 


(50) 


Haciendo un cambio a variables de estado y sus derivadas, se puede llevar 
el problema a una Ecuacion de Estado de la forma: 

X(t) = A.X(t) + B.f(t) (51) 


Donde X es el vector de estado, f el vector de input (aceleracion del suelo), 
A la matriz de estado y B la matriz de input. De esta forma en vez de resolver dos 
ecuaciones diferenciales ordinarias de segundo orden se resuelven cuatro 
ecuaciones diferenciales ordinarias de primer orden: 
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1X1 


1X1 


Xi 

X c 

d 

x c 


x c 

x t 

dt 

Xi 


fl (Xi,X c ,X v X c ) 

Xc- 


X. 


-A fait X c , X t , Xc'). 


El problems se puede resolver ahora integrando las cuatro ecuaciones 
diferenciales y definiendo cuatro condiciones iniciales de las variables buscadas. 
Asi se obtienen los dos desplazamientos (impulsivo y convectivo) del sistema y sus 
respectivas velocidades. 

2.4 Integracion numerics 

Se emplea el metodo de Runge Kutta explfcito de orden 3: 


y n +i - y n + ^ (Ai + 2 k 2 + k 3 ) 

(53) 

Kl = f(.t n ,y n ) = fn 

(54) 

h h 

(55) 

K 2 - f(t n + - , y n + -) 

K 3 = f(t n + h,y n - K t h + 2K 2 h) 

(56) 


Este metodo permite integrar numericamente las ecuaciones diferenciales. 
El siguiente valor (y n+ i) es determinado por el presente valor (y n ) mas el producto 
del tamano del intervalo (h) por una pendiente estimada. La pendiente es un 
promedio ponderado de pendientes, donde Ki es la pendiente al principio del 
intervalo; «2 es la pendiente en el punto medio del intervalo, usando Ki para 
determinar el valor de “y” en el punto t n +h/2 mediante el metodo de Euler; y K 3 es la 
pendiente al final del intervalo, con el valor de “y” determinado por K 2 . Se promedian 
las tres pendientes asignando mayor peso a la pendiente en el punto medio. 

2.5 Sistema de aislacion FPS 

Los aisladores de pendulo friccional son los mas recomendables para 
colocar en la base de tanques de almacenamiento de Ifquidos, dado que el periodo 
de aislacion Tb depende solo del radio de curvatura R c de la superficie concava de 
deslizamiento, por lo que se mantiene constante, aunque cambie el peso del tanque 
debido a variaciones en el nivel del liquido contenido. El periodo de aislacion se 
obtiene como: 


T b 



(57) 


La fuerza desarrollada en el aislador es la suma de una fuerza restitutiva, 
generada por la componente tangencial del peso W actuando sobre la superficie 
curva, y una fuerza de rozamiento, provocada por la reaccion de la componente 
normal del peso actuando sobre la superficie curva. La fuerza restitutiva controla la 
rigidez del sistema en funcion del radio de curvatura R c de la superficie concava de 
deslizamiento, y la fuerza de rozamiento controla el amortiguamiento del sistema en 
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funcion del coeficiente de friccion p entre el deslizador y la superficie concava de 
deslizamiento. 


W 

F b = — x + gW sign x (58) 

Durante el funcionamiento del aislador, la carga axial W puede variar por la 
adicion al peso de los efectos de incremento sfsmico que genera el momento de 
vuelco de la estructura, y por la componente vertical de la aceleracion del suelo. Por 
otro lado, el coeficiente de friccion p puede variar de acuerdo a la velocidad de 
deslizamiento y a la presion normal sobre la superficie de deslizamiento. En este 
artfculo, por simplicidad, no se consideran dichas variaciones y se asumen como 
valores constantes. Un analisis mas completo, parafuturas investigaciones, deberfa 
incluirlos. 

El coeficiente de friccion se incrementa rapidamente con la velocidad de 
deslizamiento hasta un cierto valor mas alia del cual se mantiene casi constante. 
Para el caso de contacto entre un deslizador de teflon y una superficie de 
deslizamiento de acero inoxidable, este efecto se explica de la siguiente forma. Para 
velocidades bajas y despues del ciclo inicial, una capa muy fina de teflon se 
transfiere al acero inoxidable y el deslizamiento ocurre ente teflon y teflon. Para altas 
velocidades, grandes escamas de teflon son removidas de la superficie y molidas 
en la interfase de deslizamiento sin transferirse al acero, lo que resulta en un 
aumento de la friccion. El aumento del coeficiente de friccion con la velocidad tiende 
a estabilizarse a velocidades del orden de 1 0 cm/s. 

Por otro lado, el coeficiente de friccion se reduce con el incremento de la 
presion de contacto sobre el deslizador hasta un valor Ifmite mas alia del cual 
permanece constante. El efecto del coeficiente de friccion con la presion es menos 
intuitivo. Al aumentar la presion entre las caras del material, el aumento del area de 
contacto es menor que el incremento en la carga normal, y por lo tanto el aumento 
en la fuerza de friccion es menor que el aumento en la carga vertical, y dado que el 
coeficiente de friccion es el cociente entre la fuerza de friccion y la carga vertical, su 
valor disminuye. 




Figura 7 Corte transversal de un aislador FPS de simple curvatura, con x=0 (izq.) 

y x=Xmax (der.) 


2.6 Analisis parametrico 

En los pasos 2.1 a 2.4 se explica el metodo utilizado para resolver el 
problema. Dichos pasos han sido programados en el software MATLAB para 
permitir una rapida aplicacion del metodo cuando se modifican los datos del 
problema. La rutina correspondiente se incluye en el Anexo A. En este artfculo se 
realiza un analisis parametrico para determinar la incidencia del perfodo de 
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aislacion, el coeficiente de friccion de los aisladores FPS y la relacion de esbeltez 
(H/R) del tanque en la respuesta sismica de este ultimo cuando tiene su base 
aislada. Tambien se analiza la efectividad del sistema de aislacion comparando la 
respuesta del tanque con base fija. Se adoptan distintos valores para los parametros 
mencionados. Para el periodo de aislacion se toman 3.0s, 3.5s, 4.0s y 4.5s. Para el 
coeficiente de friccion se toman 0.05, 0.10, 0.15 y 0.20. Para la relacion de esbeltez 
se toman 0.5, 1 .0, 1 .5 y 2.0. Se utilizan como input para el sistema cuatro registros 
sfsmicos con caracteristicas impulsivas, a fin de realizar un analisis de la historia de 
la respuesta. Los registros usados y sus parametros caracterfsticos se indican en la 
Tabla 1. 


Tabla 1 Registros sfsmicos usados en el analisis 


Lugar del 
Terremoto 

Estacion de 
registro 

Magnitud 

Momenta 

(Mw) 

PGA [g] 

PGV [m/s] 

Periodo 
pulso [s] 

Northridge, EE 
UU 

Sylmar 

Hospital 

6.7 

0.84 

1.29 

0.76 

Chichi, China 

Chichi 

7.6 

0.81 

1.26 

1.13 

Cape 

Mendocino, EE 
UU 

Petrolia 

7.0 

0.66 

0.90 

0.68 

Imperial Valley, 
EE UU 

El Centro 
array #6 

6.5 

0.44 

1.10 

1.27 


Se cuantifica la respuesta tomando como output el desplazamiento relativo 
de la masa convectiva (x c ), el desplazamiento relativo del sistema de aislacion (xb) 
que por las consideraciones realizadas en el modelo (conexion rfgida de la base del 
tanque a las paredes tambien rfgidas) es igual al desplazamiento de la masa 
impulsiva (xi), y el code basal normalizado (que puede interpretarse como el 
coeficiente sfsmico o aceleracion efectiva del sistema en la base expresada como 
fraccion de la aceleracion de la gravedad). Este ultimo parametro se obtiene para el 
caso de base aislada a partir de la expresion (58) y para el caso de base fija a partir 
de la ecuacion propuesta por Housner: 


F b = VCSai-m^ + CS^.m,) 2 (59) 

Donde S a \ y S ac son las ordenadas espectrales correspondientes al periodo 
de la masa impulsiva y convectiva, respectivamente. Por estar la masa impulsiva 
rfgidamente conectada al suelo (paredes del tanque sin flexibilidad), Sai = PGA. La 
ecuacion para Fb es conservadora, dado que la respuesta maxima de cada modo 
de vibracion (impulsivo y convectivo) no se produce en el mismo instante. 
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3 RESULTADOS 


A continuation, se presentan los valores maximos obtenidos del analisis de 
la historia de la respuesta, segun se describio en los apartados anteriores. 

3.1 Incidencia del peri'odo de aislacion en la respuesta del tanque 



Figura 8 Variacion del desplazamiento convectivo x c con el cambio de periodo de 
aislacion Tb, para un tanque ancho (H/R=0.50) - p = 0.05 
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Figura 9 Variacion del desplazamiento convectivo x c con el cambio de periodo de 
aislacion T b , para un tanque esbelto (H/R=2.00) - p = 0.05 
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Figura 10 Variation del desplazamiento del aislador Xb con el cambio de periodo 
de aislacion Tb, para un tanque ancho (H/R=0.50) - p = 0.05 
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Figura 11 Variacion del desplazamiento del aislador Xb con el cambio de periodo 
de aislacion Tb, para un tanque esbelto (H/R=2.00) - p = 0.05 
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Figura 12 Variacion del code basal normalizado Fb/W con el cambio de periodo de 
aislacion Tb, para un tanque ancho (H/R=0.50) - p = 0.05 
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Figura 13 Variacion del code basal normalizado Fb/W con el cambio de periodo de 
aislacion Tb, para un tanque esbelto (H/R=2.00) - p = 0.05 

Se observa que, en la mayoria de los registros sismicos, el aumento del 
periodo con aislacion no genera una reduccion apreciable en el desplazamiento de 
la masa convectiva del tanque, responsable del oleaje en la supedicie libre. Queda 
claramente expuesto que el aumento del periodo con aislacion conduce a una 
impodante reduccion del code basal normalizado. Por ultimo, se aprecia que, en la 
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mayorfa de los registros, a exception de Chichi, el cambio en el periodo con 
aislacion no modifica sustancialmente la respuesta en terminos de desplazamientos. 
Para el registro de Chichi se observa que, con el aumento del periodo de aislacion, 
se produce una diminution de la respuesta de desplazamientos en los tanques 
anchos y un aumento en el caso de los tanques esbeltos. 

3.2 Incidencia del coeficiente de friccion en la respuesta del tanque 
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Figura 14 Variacion del desplazamiento convectivo x c con el cambio del 
coeficiente de friccion p, para un tanque ancho (H/R=0.50) - Tb = 3.00s 
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Figura 15 Variacion del desplazamiento convectivo x c con el cambio del 
coeficiente de friccion p, para un tanque esbelto (H/R=2.00) - Tb = 3.00s 
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Figure 16 Variacion del desplazamiento del aislador Xb con el cambio del 
coeficiente de friccion p, para un tanque ancho (H/R=0.50) - Tb = 3.00s 



Figure 17 Variacion del desplazamiento del aislador Xb con el cambio del 
coeficiente de friccion p, para un tanque esbelto (H/R=2.00) - Tb = 3.00s 
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Figura 18 Variacion del code basal normalizado Fb/W con el cambio del 
coeficiente de friccion g, para un tanque ancho (H/R=0.50) - Tb = 3.00s 



Figura 19 Variacion del code basal normalizado Fb/W con el cambio del 
coeficiente de friccion p, para un tanque esbelto (H/R=2.00) - Tb = 3.00s 


Se observa que, el aumento del coeficiente de friccion provoca una 
reduccion del desplazamiento convectivo y del desplazamiento del aislador. Por otro 
lado, al incrementarse la friccion, en general aumenta el code basal normalizado, a 
excepcion del registro de Chichi en tanques anchos, y Chichi e Imperial Valley en 
tanques esbeltos, donde se observa una disminucion del code. 
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3.3 Incidencia de la relacion de esbeltez en la respuesta del tanque 
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Figura 20 Variacion del desplazamiento convectivo x c para un periodo de aislacion 

Tb=3.00s - p = 0.05 
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Figura 21 Variacion del desplazamiento convectivo x c para un periodo de aislacion 

Tb=4.50s - p = 0.05 
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Figura 22 Variacion del desplazamiento del aislador Xb para un periodo de 
aislacion Tb=3.00s - p = 0.05 
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Figura 23 Variacion del desplazamiento del aislador Xb para un periodo de 
aislacion Tb=4.50s - p = 0.05 

Se observan comportamientos diferenciados al aumentar la relacion de 
esbeltez del tanque, mientras que Chichi e Imperial Valley aumentan su respuesta, 
Northridge y Cape Mendocino la reducen. 
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3.4 Efectividad del sistema de aislacion 
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Figura 24 Variacion del code basal normalizado Fb/W con el cambio de la relacion 
de esbeltez, para un tanque no aislado y otro aislado con periodo de aislacion 
Tb=3.00s y coeficiente de friccion p=0.05 

El sistema de aislacion se muestra efectivo con todos los registros, 
tendiendo a un maximo de reduction del code basal para la relacion de esbeltez 
H/R = 1 . 


4 CONCLUSIONES 

A partir del analisis de los resultados se concluye que: 

a) La tecnica de aislacion sismica es efectiva para reducir la respuesta de 
los tanques; 
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b) El comportamiento del tanque bajo sismos impulsivos guarda estrecha 
relacion con las caracteristicas de los registros; 

c) En registros con periodos de pulso altos, como Chichi o Imperial Valley, 
la respuesta en terminos de desplazamiento se incrementa para tanques esbeltos 
en relacion a los anchos, dado que en los primeros la contribution del modo aislado 
es muy importante, y se suma a la del modo convectivo; 

d) La amplitud de potencia concentrada en bajas frecuencias en sismos de 
falla cercana afecta la respuesta de oleaje del tanque; 

e) El incremento del coeficiente de friccion reduce la efectividad del sistema 
de aislacion en terminos de aceleraciones, incrementando el code basal 
normalizado, pero ocasiona una reduccion del desplazamiento convectivo, por su 
efecto de amortiguamiento sobre el sistema; 

f) El coeficiente de friccion resulta ser un parametro clave para balancear la 
respuesta de tanques anchos, entre una reduccion de aceleracion basal y un 
aumento de desplazamiento convectivo, el valor optimo del mismo debe analizarse 
para cada registro sismico particular; 

g) El metodo propuesto permite obtener la respuesta sismica de tanques 
aislados con razonable exactitud y un muy bajo costo computacional. 
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6 ANEXO A: RUTINA EN MATLAB 

%xi desplazamiento horizontal masa impulsiva 

%xc desplazamiento horizontal masa convectiva 

%vi velocidad de la masa impulsiva 

%vc velocidad de la masa convectiva 

%ai aceleracion de la masa impulsiva 

%ac aceleracion de la masa convectiva 

%as aceleracion del suelo 

%hi altura desde la base de la masa impulsiva 

%hc altura desde la base de la masa convectiva 

%kc rigidez resorte masa convectiva 

%ec factor de amortiguamiento convectivo 

%R radio curvatura pendulo 

%u coeficiente friccion pendulo 

%g aceleracion de la gravedad 

%La distancia centro tanque a pendulo 

syms xi xc vi vc ai ac as hi he mi me kc ec R u g La A 

% Grados de libertad 
q= [xi ; xc] 
qv= [ vi ; vc] 
qa= [ai; ac] 

% Vectores de posicion de las masas 
ri= [xi ; hi+R- (R A 2-xi A 2) A 0.5] 
rc= [xi+xc; hc+R- (R A 2-xi A 2) A 0.5] 

%Vectores de velocidad 
rvi= j acobian ( ri , q) *qv 
rvc=j acobian (rc, q) *qv 

%Energia cinetica 
Ti=0 . 5*mi*rvi . ’ *rvi 
Tc=0 . 5*mc*rvc . ' *rvc 
T=Ti+Tc 

%Energia potencial gravitatoria 
Vgi=g*mi*ri (2) 

Vgc=g*mc*rc (2) 

%Energia potencial elastica 
Vsi=0 

Vsc=0 . 5*kc* (xc-xi) A 2 

%Energia potencial total 

Vi=Vgi+Vsi 

Vc=Vgc+Vsc 

V=Vi+Vc 


B real 
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%Fuerza de rozamiento pendulo en direccion tangencial 

N= ( (g*cos (asin (xi/R) ) ) + ( (1/R) * ( (videos (asin (xi/R) ) ) A 2) ) ) * ( (mi 

+mc) /2 ) 

Frt_l=-u*N*sign (vi) 

Frt_2=-u*N*sign (vi) 

%Fuerza de rozamiento pendulo en coordenadas x,y 
Fr_l = [Frt_l*cos (asin (xi/R) ) ; Frt_l*sin (asin (xi/R) ) ] 

Fr_2 = [Frt_2*cos (asin (xi/R) ) ; Frt_2*sin (as in (xi/R) ) ] 

%Vector de posicion del pendulo 
rp_l = [xi-La; R- (R A 2-xi A 2) A 0 . 5] 
rp_2= [xi+La; R- (R A 2-xi A 2 ) A 0 . 5 ] 

%Trabajo de Fr 
Fp_l=rp_l . 1 *Fr_l 
Fp~2=rp~2 . 1 *Fr~2 

%Proyeccion de Fp respecto a los grados de libertad 
Fpq=j acobian (Fp_l , q) + j acobian (Fp_2 A q) 

%Amortiguamiento convectivo 
wc= ( kc/mc) A 0 . 5 
cc=2*mc*wc*ec 


%Fuerza de amortiguamiento en coordenadas x,y 
Fa= [ -cc^vc^sign ( vc) ; 0 ] 

%Trabajo de Fa 
Fc=rc . 1 ^Fa 

%Proyeccion de Fc respecto a los grados de libertad 
Fcq=j acobian (Fc, q) 

%Fuerza no conservativa 

Qnc= [ Fpq ( 1 ) +Fcq ( 1 ) ; Fpq ( 2 ) +Fcq ( 2 ) ] 

%Primer termino ecuacion de Euler-Lagrange 
Al=j acobian (j acobian (T, qv) , qv) 

A2=j acobian (j acobian (T, q) , qv) 

A=A1 * qa+A2 * qv 

%Segundo termino ecuacion de Euler-Lagrange 
B= j acobian (T, q) 

%Tercer termino ecuacion de Euler-Lagrange 
C= j acobian (V, q) 


220 


Juan Pablo Cordone, Miguel Tornello 


%Introduccion de la fuerza sismica 
Ms= [mi, 0; 0,mc] %Matriz de masa 
n=[l;l] %Vector de colocacion 
Fs=Ms*n*as 


%Ecuacion Euler Lagrange (dos ecuaciones diferenciales de 
segundo orden) 

EL=A-B. ' +C. ' -Qnc+Fs 

%Despeje de ecuacion de segundo orden (despejo aceleracion, 
en forma manual, usar pretty (EL) para ver mejor) 

L (1,1)= (me + mi + (mc*xi A 2)/ (R A 2 - xi A 2) + (mi*xi A 2)/ (R A 2 - 
xi A 2) ) 

L ( 1 , 2 ) =mc 
L ( 2 , 1 ) =mc 
L ( 2 , 2 ) =mc 

FNL (1, 1) =EL (1) -ai*L (1, 1) -ac*L (1, 2) 

FNL (2, 1) =EL (2) -ai*L (2,1) -ac*L (2,2) 
a=-inv (L) *FNL 


%Cambio de variables (asigno variable "y" a desplazamientos y 
velocidades) : y(l)=xi y(2)=xc y(3)=vi y(4)=vc reemplazo en 
a(l,l)=ecl y en a(2,l)=ec2) 

%ecl= ( (R A 2 - y (1) A 2) * (as*mc + (kc*(2*y(2) - 2*y(l)))/2 + 
2*ec*mc*y (4) *sign (y (4) ) * (kc/mc) A (1/2) ) ) / (mi*R A 2 + mc*y(l) A 2) 

- ( (R A 2 - y (1) A 2) * (as*mi + y (3) * ( (2*mc*y (3) *y (1) ) / (R A 2 - 
y(l) A 2) + (2*mi*y(3) *y(l) ) / (R A 2 - y(l) A 2) + 

(2*mc*y (3) *y (1) A 3) / (R A 2 - y(l) A 2) A 2 + (2*mi*y (3) *y (1) A 3) / (R A 2 

- y(l) A 2) A 2) - (kc* (2*y (2) - 2*y(l)))/2 - 

(mc*y (3) A 2*y (1) A 3) / (R A 2 - y(l) A 2) A 2 - (mi*y (3) A 2*y (1) A 3) / (R A 2 

- y ( 1 ) A 2 ) A 2 + (g*mc*y (1) ) / (R A 2 - y (1) A 2) A (1/2) + 

(g*mi*y (1) ) / (R A 2 - y (1) A 2) A (1/2) - (mc*y (3) A 2*y (1) ) / (R A 2 - 
y ( 1 ) A 2 ) - (mi*y (3) A 2*y (1) ) / (R A 2 - y(l) A 2) + 2*u*sign (y (3) ) * (1 

- y (1) A 2/R A 2) A (1/2) * (mc/2 + mi/2)*(g*(l - y (1) A 2/R A 2) A (1/2) - 
(y (3) A 2* (y (1) A 2/R A 2 - 1))/R) + 

2*ec*mc*y (4) *sign (y (4) ) * (kc/mc) A (1/2) + (2*u*sign (y (3) ) * (R - 
(R A 2 - y (1) A 2) A (1/2) )* (mc/2 + mi/2)*(g*(l - y (1) A 2/R A 2) A (1/2) 

- (y (3) A 2* (y (1) A 2/R A 2 - 1))/R))/R + u*sign (y (3) ) * (1 - 

y(l) A 2/R A 2) A (1/2) * ( (2*y(3) A 2*y(l) ) /R A 3 + (g*y (1) ) / (R A 2* (1 - 
y (1) A 2/R A 2) A (1/2) ) ) * (La - y(l))*(mc/2 + mi/2) - 
u*sign (y (3) ) * (La + y(l))*(l - 

y (1) A 2/R A 2) A (1/2) * ( (2*y (3) A 2*y (1) ) /R A 3 + (g*y (1) ) / (R A 2* (1 - 
y (1) A 2/R A 2) A (1/2) ))* (mc/2 +mi/2) + 

(2*u*y (1) A 2*sign(y(3) ) * (mc/2 + mi/2)*(g*(l - 
y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y (1) A 2/R A 2 - 1 ) ) /R) ) / (R* (R A 2 - 
y (1) A 2) A (1/2) ) - (2*u*y(l)*sign(y(3))*((2*y(3) A 2*y(l))/R A 3 + 
(g*y (1) ) / (R A 2* (1 - y (1) A 2/R A 2) A (1/2) ) ) * (R - (R A 2 - 
y (1) A 2) A (1/2) )* (mc/2 + mi/2))/R + (u*y (1) *sign (y (3) ) * (La - 
y(l))*(mc/2 + mi/2 ) * (g* ( 1 - y (1) A 2/R A 2) A (1/2) - 
(y (3) A 2* (y (1) A 2/R A 2 - 1) ) /R) ) / (R A 2* (1 - y (1) A 2/R A 2) A (1/2) ) - 

(u*y (1) *sign (y (3) ) * (La + y(l))*(mc/2 + mi/2)*(g*(l - 
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y (1) A 2/R A 2) * (1/2) - (y(3) A 2* (y (1) A 2/R A 2 - 1) ) /R) ) / (R A 2* (1 - 
y (1) A 2/R A 2) A (1/2) ) ) ) / (mi*R A 2 + mc*y(l) A 2) 

%ec2= ( (R A 2 - y (1) A 2) * (as*mi + y (3) * ( (2*mc*y (3) *y (1) ) / (R A 2 - 
y ( 1 ) A 2 ) + (2*mi*y (3) *y (1) ) / (R A 2 - y(l) A 2) + 

(2*mc*y (3) *y (1) A 3) / (R A 2 - y(l) A 2) A 2 + (2*mi*y (3) *y (1) A 3) / (R A 2 

- y(l) A 2) A 2) - (kc* (2*y (2) - 2*y(l)))/2 - 

(mc*y (3) A 2*y (1) A 3) / (R A 2 - y(l) A 2) A 2 - (mi*y (3) A 2*y (1) A 3) / (R A 2 

- y(l) A 2) A 2 + (g*mc*y(l) ) / (R A 2 - y (1) A 2) A (1/2) + 

(g*mi*y (1) ) / (R A 2 - y (1) A 2) A (1/2) - (mc*y (3) A 2*y (1) ) / (R A 2 - 
y ( 1 ) A 2 ) - (mi*y (3) A 2*y (1) ) / (R A 2 - y(l) A 2) + 2*u*sign (y (3) ) * (1 

- y (1) A 2/R A 2) A (1/2) * (mc/2 + mi/2)*(g*(l - y (1) A 2/R A 2) A (1/2) - 
(y (3) A 2* (y (1) A 2/R A 2 - 1))/R) + 

2*ec*mc*y (4) *sign (y (4) ) * (kc/mc) A (1/2) + (2*u*sign (y (3) ) * (R - 
(R A 2 - y (1) A 2) A (1/2) )* (mc/2 + mi/2 ) * (g* (1 - y (1) A 2/R A 2) A (1/2) 

- (y (3) A 2* (y (1) A 2/R A 2 - 1))/R))/R + u*sign (y (3) ) * (1 - 

y (1) A 2/R A 2) A (1/2) * ( (2*y (3) A 2*y (1) ) /R A 3 + (g*y (1) ) / (R A 2* (1 - 
y (1) A 2/R A 2) A (1/2) ) ) * (La - y(l))*(mc/2 + mi/2) - 
u*sign (y (3) ) * (La + y(l))*(l - 

y(l) A 2/R A 2) A (1/2) * ( (2*y(3) A 2*y(l) ) /R A 3 + (g*y (1) ) / (R A 2* (1 - 
y (1) A 2/R A 2) A (1/2) ))* (mc/2 +mi/2) + 

(2*u*y (1) A 2*sign (y (3) ) * (mc/2 + mi/2)*(g*(l - 

y (1) A 2/R A 2) A (1/2) - (y(3) A 2* (y (1) A 2/R A 2 - 1 ) ) /R) ) / (R* (R A 2 - 

y (1) A 2) A (1/2) ) - (2*u*y(l) *sign(y(3) ) * ( (2*y(3) A 2*y(l) ) /R A 3 + 

(g*y (1) ) / (R A 2* (1 - y (1) A 2/R A 2) A (1/2) ) ) * (R - (R A 2 - 
y (1) A 2) A (1/2) )* (mc/2 + mi/2))/R + (u*y (1) *sign (y (3) ) * (La - 
y(l))*(mc/2 + mi/ 2 ) * (g* ( 1 - y (1) A 2/R A 2) A (1/2) - 
(y (3) A 2* (y (1) A 2/R A 2 - 1) ) /R) ) / (R A 2* (1 - y (1) A 2/R A 2) A (1/2) ) - 

(u*y (1) *sign (y (3) ) * (La + y(l))*(mc/2 + mi/2)*(g*(l - 
y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y (1) A 2/R A 2 - 1) ) /R) ) / (R A 2* (1 - 

y (1) A 2/R A 2) A (1/2) ) ) ) / (mi*R A 2 + mc*y(l) A 2) - ( (R A 2*mc + 

R A 2*mi) * (as*mc + (kc*(2*y(2) - 2*y(l)))/2 + 

2*ec*mc*y (4) *sign(y(4) ) * (kc/mc) A (1/2) ) ) / (me* (mi*R A 2 + 
mc*y ( 1 ) A 2 ) ) 


%Creacion de funcion para integrar numericamente 
f=@ (t,y,as) [y (3) ;y(4) ; ( (R A 2 - y ( 1 ) A 2 ) * (as*mc + (kc*(2*y(2) - 
2*y (1) ) ) / 2 + 2*ec*mc*y(4) *sign(y(4) )* (kc/mc) A (1/2) ))/ (mi*R A 2 
+ mc*y ( 1 ) A 2 ) - ( (R A 2 - y ( 1 ) A 2 ) * (as*mi + 

y (3) * ( (2*mc*y (3) *y (1) ) / (R A 2 - y(l) A 2) + (2*mi*y (3) *y (1) ) / (R A 2 
- y ( 1 ) A 2 ) + (2*mc*y (3) *y (1) A 3) / (R A 2 - y(l) A 2) A 2 + 

(2*mi*y (3) *y (1) A 3) / (R A 2 - y(l) A 2) A 2) - (kc*(2*y(2) - 

2*y (1) ) ) / 2 - (mc*y (3) A 2*y (1) A 3) / (R A 2 - y(l) A 2) A 2 - 
(mi*y (3) A 2*y (1) A 3) / (R A 2 - y(l) A 2) A 2 + (g*mc*y ( 1 ) ) / (R A 2 - 
y (1) A 2) A (1/2) + (g*mi*y (1) ) / (R A 2 - y (1) A 2) A (1/2) - 
(mc*y (3) A 2*y (1) ) / (R A 2 - y(l) A 2) - (mi*y (3) A 2*y (1) ) / (R A 2 - 
y(l) A 2) + 2*u*sign (y (3) ) * (1 - y (1) A 2/R A 2) A (1/2) * (mc/2 + 
mi/ 2 ) * (g* ( 1 - y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y ( 1 ) A 2/R A 2 - 
1))/R) + 2*ec*mc*y (4) *sign (y (4) ) * (kc/mc) A (1/2) + 

(2*u*sign (y (3) ) * (R - (R A 2 - y (1) A 2) A (1/2) )* (mc/2 + 
mi/ 2 ) * (g* ( 1 - y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y (1) A 2/R A 2 - 

1))/R))/R + u*sign (y (3) ) * (1 - 

y (1) A 2/R A 2) A (1/2) * ( (2*y (3) A 2*y (1) ) /R A 3 + (g*y (1) ) / (R A 2* (1 - 
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y (1) A 2/R A 2) A (1/2) ) ) * (La - y(l))*(mc/2 + mi/2) - 
u*sign (y (3) ) * (La + y(l))*(l - 

y(l) A 2/R A 2) a (1/2) * ( (2*y(3) A 2*y(l) ) /R A 3 + (g*y (1) ) / (R A 2* (1 - 
y (1) A 2/R A 2) A (1/2) ) ) * (mc/2 +mi/2) + 

(2*u*y(l) A 2*sign(y(3) ) * (mc/2 + mi/2)*(g*(l - 

y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y (1) A 2/R A 2 - 1 ) ) /R) ) / (R* (R A 2 - 

y (1) A 2) A (1/2) ) - (2*u*y(l) *sign (y(3) ) * ( (2*y(3) A 2*y(l) ) /R A 3 + 

(g*y (1) ) / (R A 2* (1 - y (1) A 2/R A 2) A (1/2) ) ) * (R - (R A 2 - 

y (1) A 2) A (1/2) )* (mc/2 + mi/2))/R + (u*y (1) *sign (y (3) ) * (La - 

y ( 1 ) ) * (me/ 2 + mi/ 2 ) * (g* ( 1 - y (1) A 2/R A 2) A (1/2) - 

(y (3) A 2* (y (1) A 2/R A 2 - 1) ) /R) ) / (R A 2* (1 - y ( 1 ) A 2/R A 2 ) A ( 1/2 ) ) - 

(u*y (1) *sign (y (3) ) * (La + y(l))*(mc/2 + mi/2)*(g*(l - 

y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y (1) A 2/R A 2 - 1) ) /R) ) / (R A 2* (1 - 

y (1) A 2/R A 2) A (1/2) ) ) ) / (mi*R A 2 +mc*y(l) A 2); ( (R A 2 - 

y (1) A 2) * (as*mi + y (3) * ( (2*mc*y (3) *y (1) ) / (R A 2 - y(l) A 2) + 

(2*mi*y (3) *y (1) ) / (R A 2 - y(l) A 2) + (2*mc*y (3) *y ( 1 ) A 3) / (R A 2 - 

y(l) A 2) A 2 + (2*mi*y(3) *y(l) A 3) / (R A 2 - y(l) A 2) A 2) - 

(kc* (2*y (2) - 2*y (1) ) ) / 2 - (mc*y (3) A 2*y ( 1 ) A 3) / (R A 2 - 

y ( 1 ) A 2 ) A 2 - (mi*y(3) A 2*y(l) A 3) / (R A 2 - y(l) A 2) A 2 + 

(g*mc*y (1) ) / (R A 2 - y (1) A 2) A (1/2) + (g*mi*y ( 1 ) ) / (R A 2 - 
y (1) A 2) A (1/2) - (mc*y(3) A 2*y(l) ) / (R A 2 - y(l) A 2) - 

(mi*y (3) A 2*y (1) ) / (R A 2 - y(l) A 2) + 2*u*sign (y (3) ) * (1 - 
y(l) A 2/R A 2) A (1/2) * (mc/2 + mi/2)*(g*(l - y ( 1 ) A 2/R A 2 ) A ( 1/2 ) - 
(y (3) A 2* (y (1) A 2/R A 2 - 1))/R) + 

2*ec*mc*y (4) *sign (y (4) ) * (kc/mc) A (1/2) + (2*u*sign (y (3) ) * (R - 
(R A 2 - y (1) A 2) A (1/2) )* (mc/2 + mi/2)*(g*(l - y ( 1 ) A 2/R A 2 ) A ( 1/2 ) 
- (y (3) A 2* (y (1) A 2/R A 2 - 1))/R))/R + u*sign (y (3) ) * (1 - 
y (1) A 2/R A 2) A (1/2) * ( (2*y (3) A 2*y (1) ) /R A 3 + (g*y (1) ) / (R A 2* (1 - 
y (1) A 2/R A 2) A (1/2) ) ) * (La - y(l))*(mc/2 + mi/2) - 
u*sign (y (3) ) * (La + y(l))*(l - 

y(l) A 2/R A 2) A (1/2) * ( (2*y(3) A 2*y(l) ) /R A 3 + (g*y (1) ) / (R A 2* (1 - 
y (1) A 2/R A 2) A (1/2) ))* (mc/2 +mi/2) + 

(2*u*y (1) A 2*sign (y (3) ) * (mc/2 + mi/2)*(g*(l - 

y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y (1) A 2/R A 2 - 1 ) ) /R) ) / (R* (R A 2 - 

y (1) A 2) A (1/2) ) - (2*u*y(l) *sign (y(3) ) * ( (2*y(3) A 2*y(l) ) /R A 3 + 

(g*y (1) ) / (R A 2* (1 - y (1) A 2/R A 2) A (1/2) ) ) * (R - (R A 2 - 
y (1) A 2) A (1/2) )* (mc/2 + mi/2))/R + (u*y (1) *sign (y (3) ) * (La - 
y ( 1 ) ) * (me/ 2 + mi/ 2 ) * (g* ( 1 - y (1) A 2/R A 2) A (1/2) - 
(y (3) A 2* (y (1) A 2/R A 2 - 1) ) /R) ) / (R A 2* (1 - y ( 1 ) A 2/R A 2 ) A ( 1/2 ) ) - 
(u*y (1) *sign (y (3) ) * (La + y(l))*(mc/2 + mi/2)*(g*(l - 
y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y (1) A 2/R A 2 - 1) ) /R) ) / (R A 2* (1 - 

y (1) A 2/R A 2) A (1/2) ) ) ) / (mi*R A 2 + mc*y(l) A 2) - ( (R A 2*mc + 

R A 2*mi) * (as*mc + (kc*(2*y(2) - 2*y(l)))/2 + 

2*ec*mc*y (4) *sign(y(4) ) * (kc/mc) A (1/2) ) ) / (me* (mi*R A 2 + 
mc*y ( 1 ) A 2 ) ) ] 


%Asignacion de valores 

D=2 %Diametro del tanque (dato) 

H=2 %Altura de agua del tanque (dato) 
p=1000 %Densidad del liquido contenido (dato) 
g=9.81 %Aceleracion de la gravedad 
M=p*pi*H*D*D/4 %Masa total del tanque 
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hi=0.375*H %altura de la masa impulsiva 

hc=H* (1- (D/ (3 . 67*H) ) *tanh (1 . 835*H/D) ) %altura de la masa 
convectiva 

mi=M*tanh (0 . 866*D/H) / (0 . 866*D/H) %masa impulsiva 
mc=0 . 835*M*tanh (3 . 67*H/D) / (3 . 67*H/D) %masa convectiva 
kc=16 . 13042*mc*mc*g*H/ (M*D*D) %rigidez convectiva 
ec=0.005 %factor de amortiguamiento convectivo (dato) 
u=0.05 %coeficiente friccion pendulo (dato) 

R=5.03 %radio curvatura pendulo (dato) 5.03 3.98 3.04 2.24 
La=l %distancia centro tanque a pendulo 
Regis tro=load ( 'CHICHI-0 . 005-g. txt ' ) ; 

as= [Registro ( : , 2 ) * 9 . 81 ; zeros ( 18000 , 1 )]; %Aceleracion del 
suelo 

t=Registro ( : , 1) ; 

%Integracion RK3 

f=@ (t,y,as) [y (3) ;y(4) ; ( (R A 2 - y ( 1 ) A 2 ) * (as*mc + (kc*(2*y(2) - 
2*y(l)))/2 + 2*ec*mc*y(4) *sign(y(4) ) * (kc/mc) A (1/2) ) ) / (mi*R A 2 
+ mc*y ( 1 ) A 2 ) - ( (R A 2 - y ( 1 ) A 2 ) * (as*mi + 

y (3) * ( (2*mc*y (3) *y (1) ) / (R A 2 - y(l) A 2) + (2*mi*y (3) *y (1) ) / (R A 2 
- y ( 1 ) A 2 ) + (2*mc*y (3) *y (1) A 3) / (R A 2 - y(l) A 2) A 2 + 

(2*mi*y (3) *y (1) A 3) / (R A 2 - y(l) A 2) A 2) - (kc*(2*y(2) - 

2*y (1) ) ) / 2 - (mc*y (3) A 2*y (1) A 3) / (R A 2 - y(l) A 2) A 2 - 
(mi*y (3) A 2*y (1) A 3) / (R A 2 - y(l) A 2) A 2 + (g*mc*y ( 1 ) ) / (R A 2 - 
y (1) A 2) A (1/2) + (g*mi*y(l) ) / (R A 2 - y (1) A 2) A (1/2) - 

(mc*y (3) A 2*y (1) ) / (R A 2 - y(l) A 2) - (mi*y (3) A 2*y ( 1 ) ) / (R A 2 - 
y(l) A 2) + 2*u*sign(y(3) ) * (1 - y (1) A 2/R A 2) A (1/2) * (mc/2 + 
mi/ 2 ) * (g* ( 1 - y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y ( 1 ) A 2/R A 2 - 
1))/R) + 2*ec*mc*y(4) *sign(y(4) )* (kc/mc) A (1/2) + 

(2*u*sign (y (3) ) * (R - (R A 2 - y (1) A 2) A (1/2) )* (mc/2 + 
mi/ 2 ) * (g* ( 1 - y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y (1) A 2/R A 2 - 

1))/R))/R + u*sign (y (3) ) * (1 - 

y (1) A 2/R A 2) A (1/2) * ( (2*y (3) A 2*y (1) ) /R A 3 + (g*y (1) ) / (R A 2* (1 - 
y (1) A 2/R A 2) A (1/2) ) ) * (La - y(l))*(mc/2 + mi/2) - 
u*sign (y (3) ) * (La + y(l))*(l - 

y(l) A 2/R A 2) A (1/2) * ( (2*y(3) A 2*y(l) ) /R A 3 + (g*y (1) ) / (R A 2* (1 - 
y (1) A 2/R A 2) A (1/2) ))* (mc/2 +mi/2) + 

(2*u*y(l) A 2*sign(y(3) ) * (mc/2 + mi/2)*(g*(l - 

y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y (1) A 2/R A 2 - 1 ) ) /R) ) / (R* (R A 2 - 

y (1) A 2) A (1/2) ) - (2*u*y(l) *sign (y(3) ) * ( (2*y(3) A 2*y(l) ) /R A 3 + 

(g*y (1) ) / (R A 2* (1 - y (1) A 2/R A 2) A (1/2) ) ) * (R - (R A 2 - 

y (1) A 2) A (1/2) )* (mc/2 + mi/2))/R + (u*y (1) *sign (y (3) ) * (La - 

y ( 1 ) ) * (me/ 2 + mi/ 2 ) * (g* ( 1 - y (1) A 2/R A 2) A (1/2) - 

(y (3) A 2* (y (1) A 2/R A 2 - 1) ) /R) ) / (R A 2* (1 - y ( 1 ) A 2/R A 2 ) A ( 1/2 ) ) - 

(u*y (1) *sign (y (3) ) * (La + y(l))*(mc/2 + mi/2)*(g*(l - 

y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y (1) A 2/R A 2 - 1) ) /R) ) / (R A 2* (1 - 

y (1) A 2/R A 2) A (1/2) ) ) ) / (mi*R A 2 +mc*y(l) A 2); ( (R A 2 - 

y (1) A 2) * (as*mi + y (3) * ( (2*mc*y (3) *y (1) ) / (R A 2 - y(l) A 2) + 

(2*mi*y (3) *y (1) ) / (R A 2 - y(l) A 2) + (2*mc*y (3) *y ( 1 ) A 3) / (R A 2 - 

y(l) A 2) A 2 + (2*mi*y(3) *y(l) A 3) / (R A 2 - y(l) A 2) A 2) - 

(kc* (2*y (2) - 2*y (1) ) ) / 2 - (mc*y (3) A 2*y ( 1 ) A 3) / (R A 2 - 

y ( 1 ) A 2 ) A 2 - (mi*y (3) A 2*y (1) A 3) / (R A 2 - y(l) A 2) A 2 + 
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(g*mc*y (1) ) / (R A 2 - y (1) A 2) A (1/2) + (g*mi*y ( 1 ) ) / (R A 2 - 
y (1) A 2) A (1/2) - (mc*y(3) A 2*y (1) ) / (R A 2 - y(l) A 2) - 

(mi*y (3) A 2*y (1) ) / (R A 2 - y(l) A 2) + 2*u*sign (y (3) ) * (1 - 
y (1) A 2/R A 2) A (1/2) * (mc/2 + mi/2 ) * (g* (1 - y ( 1 ) A 2/R A 2 ) A ( 1/2 ) - 
(y (3) A 2* (y (1) A 2/R A 2 - 1))/R) + 

2*ec*mc*y (4) *sign (y (4) ) * (kc/mc) A (1/2) + (2*u*sign (y (3) ) * (R - 
(R A 2 - y (1) A 2) A (1/2) )* (mc/2 + mi/2 ) * (g* (1 - y (1) A 2/R A 2) A (1/2) 
- (y (3) A 2* (y (1) A 2/R A 2 - 1))/R))/R + u*sign (y (3) ) * (1 - 
y (1) A 2/R A 2) A (1/2) * ( (2*y (3) A 2*y (1) ) /R A 3 + (g*y (1) ) / (R A 2* (1 - 
y (1) A 2/R A 2) A (1/2) ) ) * (La - y(l))*(mc/2 + mi/2) - 
u*sign (y (3) ) * (La + y(l))*(l - 

y (1) A 2/R A 2) A (1/2) * ( (2*y (3) A 2*y (1) ) /R A 3 + (g*y (1) ) / (R A 2* (1 - 
y (1) A 2/R A 2) A (1/2) ))* (mc/2 +mi/2) + 

(2*u*y (1) A 2*sign (y (3) ) * (mc/2 + mi/2)*(g*(l - 

y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y (1) A 2/R A 2 - 1 ) ) /R) ) / (R* (R A 2 - 

y (1) A 2) A (1/2) ) - (2*u*y(l) *sign (y (3) ) * ( (2*y(3) A 2*y(l) ) /R A 3 + 

(g*y (1) ) / (R A 2* (1 - y (1) A 2/R A 2) A (1/2) ) ) * (R - (R A 2 - 
y (1) A 2) A (1/2) )* (mc/2 + mi/2))/R + (u*y (1) *sign (y (3) ) * (La - 
y ( 1 ) ) * (me/ 2 + mi/ 2 ) * (g* ( 1 - y (1) A 2/R A 2) A (1/2) - 
(y (3) A 2* (y (1) A 2/R A 2 - 1) ) /R) ) / (R A 2* (1 - y (1) A 2/R A 2) A (1/2) ) - 

(u*y (1) *sign (y (3) ) * (La + y(l))*(mc/2 + mi/2)*(g*(l - 
y (1) A 2/R A 2) A (1/2) - (y (3) A 2* (y (1) A 2/R A 2 - 1) ) /R) ) / (R A 2* (1 - 

y (1) A 2/R A 2) A (1/2) ) ) ) / (mi*R A 2 + mc*y(l) A 2) - ( (R A 2*mc + 

R A 2*mi) * (as*mc + (kc*(2*y(2) - 2*y(l)))/2 + 

2*ec*mc*y (4) *sign(y(4) ) * (kc/mc) A (1/2) ) ) / (me* (mi*R A 2 + 
mc*y(l) A 2 ) ) ] 

tO = 0; % Ingresar condicion inicial de t 

yO = [0;0;0;0] ; % Ingresar condicion inicial de y 

T = 150; % Ingresar valor final de t 

h = 0.005; % Ingresar tamano de paso de integracion 

N = round ( (T-t0)/h ); % Numero de pasos desde tO a T; 

ROUND (X) redondea los elementos de X al entero mas proximo 
tn = linspace (tO, T, N+l) ; % Intervalo de tiempo; 

LINSPACE ( tO , T, N+l) genera N+l puntos entre tO y T 
y = zeros (4, N+l); % Intervalo de vector de estado; ZEROS (1, 
N+l) es una matriz de ceros de 4 filas y (N+l) columnas 
y(:, 1) = yO; % Almacenamiento del valor incial xO en x 
t = tO ; % Inicializa t con tO 

for i = 1:N % Ciclo de Runge-Kutta (explicito de 1 paso, 
necesita conocer 1 solo punto anterior) para hallar los 2 
primeros puntos dato necesarios para Adams-Bashf orth 
(explicito de multiple paso); explicito: usa informacion 
hacia atras (conocida) 

kl = f (t, y ( : , i) , as (i+1) ) ; % Evaluacion del incremento de 
pendiente al comienzo del intervalo usando y' (Euler) 

k2 = f(t+h/2, y ( : , i) +h/2*kl, as (i+1)); % Evaluacion del 
incremento de pendiente a mitad del intervalo usando 
y'+h/2*kl 

k3 = f (t+h, y ( : , i) -h*kl+2*h*k2, as (i+1)); % Evaluacion 
del incremento de pendiente al final del intervalo usando y'- 
h*kl+2*h*k2 
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y (:,i + l) = y ( : , i) +h/4* (kl + 2*k2 + k3) ; % Determina el 
valor de y, usando el valor anterior (1 paso) mas el promedio 
ponderado de tres incrementos, donde cada incremento es el 
producto del tamano del intervalo, h , y una pendiente 
estimada especificada por la funcion f en el lado derecho de 
la ecuacion diferencial, una vez obtenido el valor de y, se 
descartan las 3 evaluaciones realizadas 

t = t + h; % Determina el punto de tiempo a evaluar en el 
siguiente ciclo 
end 

desplimax=max (abs (y ( 1 , : ) ) ) % Desplazamiento impulsivo maximo 

desplcmax=max (abs (y (2 , : ) ) ) % Desplazamiento convectivo maximo 

Cs=max (abs (y ( 1 , : ) . /R+sign (y ( 3 , : ) ) . *u) ) % Corte basal 

normalizado maximo 


ttran=tn ' ; 
ytran=y ’ ; 
figure; 

plot (tn, y(l,:), 1 r- 1 ) ; % Grafica desplazamiento impulsivo 

xlabel ( ' t (s) ’ ) ; 
ylabel ( ' xi (m) 1 ) ; 
hold on; 

plot(tn A y(2,:), 'b-'); % Grafica desplazamiento convectivo 

xlabel ( 1 1 (s) ' ) ; 
ylabel ( ' xc (m) ’ ) ; 


